---
title: "Live Imaging PD Project"
output:   
  html_document:
   code_folding: hide
---

```{r toolbox, include=FALSE}
# Load packages

# library(gt)
# library(tidyverse)
library(ggplot2)
library(MASS)
library(lme4)
library(viridis)
library(lmtest)
library(cowplot)
# library(pscl)
# library(boot)
library(moments)
# library(lmtest)
library(dplyr) 
library(ggpubr)
# library(GLMMadaptive)
# library(MASS)
library(knitr)
# library(effectsize)
library(ggsignif)
library(lmtest)
library(glmmTMB)
library(bbmle)
library(emmeans)
```

```{r import}
#Import CSV files
live_dat.pre <- read.csv("/Users/adambindas/Dropbox/Parkinson's disease ENS Paper/Data (Raw and Organized)/Live Imaging ENS PFF/400 ms live imaging.csv", header=TRUE, stringsAsFactors=FALSE)

#Subset out replicate 19 (Butyrate dosage was altered)
live_dat <- subset(live_dat.pre, rat.no != 19)

#Factor and variable development
live_dat$Groups.f <- factor(live_dat$group, levels = c("Control","But", "PFF", "PFF+But","2XPFF", "4XPFF", "4XPFF+But", "6XPFF"))
live_dat$control <- ifelse(live_dat$group== "Control", 1, 0)
live_dat$but <- ifelse(live_dat$group== "But", 1, 0)
live_dat$pff <- ifelse(live_dat$group== "PFF", 1, 0)
live_dat$pff.but <- ifelse(live_dat$group== "PFF+But", 1, 0)
live_dat$pff.2 <- ifelse(live_dat$group== "2XPFF", 1, 0)
live_dat$pff.4 <- ifelse(live_dat$group== "4XPFF", 1, 0)
live_dat$pff.4.but <- ifelse(live_dat$group== "4XPFF+But", 1, 0)
live_dat$pff.6 <- ifelse(live_dat$group== "6XPFF", 1, 0)
live_dat$day.no <- factor(live_dat$day)
live_dat$rat.no <- factor(live_dat$rat.no)
live_dat$w.id <- paste(live_dat$rat.no,"-",live_dat$well)
live_dat$present <- ifelse(live_dat$avg.int < 1.05,1,0)
live_dat$c.avgint <- ifelse(live_dat$avg.int-1.05 < 0, 0, live_dat$avg.int-1.05)

#Final analysis subsetting
live_dat21 <- subset(live_dat,day == 21)
live_dat.21.pff <- subset(live_dat21, 
                     (Groups.f == "Control") | 
                     (Groups.f == "PFF") | 
                     (Groups.f == "2XPFF") | 
                     (Groups.f == "4XPFF") | 
                     (Groups.f == "6XPFF"))
live_dat.t <- subset(live_dat, 
                     (Groups.f == "Control") | 
                     (Groups.f == "PFF") | 
                     (Groups.f == "But") | 
                     (Groups.f == "PFF+But"))
live_dat2 <- subset(live_dat, control == 0 & but == 0 &  rat.no != 19)
live_dat2simp <- subset(live_dat, 
                          day.no == 21 &
                          (Groups.f == "Control" | Groups.f == "PFF" | Groups.f == "But"))
live_dat.t.pff <- subset(live_dat.t, Groups.f == "PFF")

# Identify sample size that is analyzed. N = experimental replicate, m = wells, d = number of images
# length(unique(g_cone$rat))

# For each level of information, identify the number of unique values (e.g., attempt number) across each day and experimental group
n.time <- aggregate(data = live_dat.t,
                          rat.no ~ Groups.f + day,
                          function(rat.no) length(unique(rat.no)))
m.time <- aggregate(data = live_dat.t,
                          w.id ~ Groups.f + day,
                          function(w.id) length(unique(w.id)))
i.time  <- aggregate(data = live_dat.t,
                          image.no ~ Groups.f + day,
                          function(image.no) length(unique(image.no)))

# Merge output sample size variables into a table form
live.time.replicates <- cbind(n.time[,2],n.time[,3],m.time[,3],i.time[,3])

# Label column names
colnames(live.time.replicates) <- c("Day","N","m","no.images")

# Label row names with a combined variable
rownames(live.time.replicates) <- m.time[,1]

# Display sample size
live.time.replicates

```

```{r day 21 combined plot}
# Identify sample size that is analyzed. N = experimental replicate, m = wells, d = number of images
# length(unique(g_cone$rat))

# For each level of information, identify the number of unique values (e.g., attempt number) across each day and experimental group
n.d21 <- aggregate(data = live_dat21,
                          rat.no ~ Groups.f,
                          function(rat.no) length(unique(rat.no)))
m.d21 <- aggregate(data = live_dat21,
                          w.id ~ Groups.f,
                          function(w.id) length(unique(w.id)))
i.d21  <- aggregate(data = live_dat21,
                          image.no ~ Groups.f,
                          function(image.no) length(unique(image.no)))

# Merge output sample size variables into a table form
live.replicates <- cbind(n.d21[,2],m.d21[,2],i.d21[,2])

# Label column names
colnames(live.replicates) <- c("N","m","no.images")

# Label row names with a combined variable
rownames(live.replicates) <-m.d21[,1]

# Display sample size
live.replicates

# Plot values
plot0 <- ggplot(live_dat21, aes(x=Groups.f, y=avg.int, fill=Groups.f)) +
  ylab("Average intensity") + xlab("Groups") + 
  labs(title = "Neural cultures by 21 respond to butyrate dosage", subtitle = "400 ms exposure time", fill = "Groups")+
  geom_bar(position = 'dodge', stat = 'summary', fun = 'mean') +
    stat_summary(fun.data = mean_se,  
                 geom = "errorbar", position = 'dodge', colour = "red") + 
  geom_signif(comparisons=list(c("Control", "But")), annotations="NS",
              y_position = 3, tip_length = 0.05, vjust=0) +
  geom_signif(comparisons=list(c("Control", "PFF")), annotations="***",
              y_position = 8, tip_length = 0.05, vjust=0.4) +
  geom_signif(comparisons=list(c("But", "PFF")), annotations="***",
              y_position = 6.5, tip_length = 0.05, vjust=0.4) +
  geom_signif(comparisons=list(c("PFF", "PFF+But")), annotations="***",
              y_position = 6.5, tip_length = 0.05, vjust=0.4) +
  geom_signif(comparisons=list(c("PFF", "2XPFF")), annotations="**",
              y_position = 8, tip_length = 0.05, vjust=0.4) +
  geom_signif(comparisons=list(c("4XPFF", "4XPFF+But")), annotations="***",
              y_position = 10.5, tip_length = 0.05, vjust=0.4) +
  geom_signif(comparisons=list(c("2XPFF", "4XPFF")), annotations="**",
              y_position = 10.5, tip_length = 0.05, vjust=0.4) +
  geom_signif(comparisons=list(c("4XPFF", "6XPFF")), annotations="*",
              y_position = 14, tip_length = 0.05, vjust=0) +
  geom_dotplot(data = live_dat21, aes(x = Groups.f, y = avg.int),
      alpha = .8, position_dodge(width = 0.9), binaxis = "y", dotsize = 0.3,
      stackdir = "center", stackratio = 1) +
  theme(legend.position="none")

plot0

plot1 <- ggplot(live_dat2, aes(x=Groups.f, y=avg.int, fill=Groups.f)) +
  ylab("Average intensity") + xlab("Groups") + 
  labs(title = "Neural cultures by 21 respond to butyrate dosage", subtitle = "400 ms exposure time", fill = "Groups")+
  geom_bar(position = 'dodge', stat = 'summary', fun = 'mean') +
    stat_summary(fun.data = mean_se,  
                 geom = "errorbar", position = 'dodge', colour = "red") + 
  geom_signif(comparisons=list(c("PFF", "PFF+But")), annotations="*",
              y_position = 6.5, tip_length = 0.05, vjust=0.4) +
  geom_signif(comparisons=list(c("PFF", "2XPFF")), annotations="*",
              y_position = 8, tip_length = 0.05, vjust=0.4) +
  geom_signif(comparisons=list(c("4XPFF", "4XPFF+But")), annotations="**",
              y_position = 10.5, tip_length = 0.05, vjust=0.4) +
  geom_signif(comparisons=list(c("2XPFF", "4XPFF")), annotations="**",
              y_position = 10.5, tip_length = 0.05, vjust=0.4) +
  geom_signif(comparisons=list(c("4XPFF", "6XPFF")), annotations="NS",
              y_position = 17.5, tip_length = 0.05, vjust=0) +
  geom_dotplot(data = live_dat2, aes(x = Groups.f, y = avg.int),
      alpha = .8, position_dodge(width = 0.9), binaxis = "y", dotsize = 0.6,
      stackdir = "center", stackratio = 1) +
  theme(legend.position="none")

plot1

live_dat2simp <- subset(live_dat, 
                          day.no == 21 &
                          (Groups.f == "Control" | Groups.f == "PFF" | Groups.f == "But"))

plot2 <- ggplot(live_dat2simp, aes(x=Groups.f, y=avg.int, fill=Groups.f)) +
  ylab("Average intensity") + xlab("Groups") + 
  labs(title = "Preformed fibril fluorescent tag is measurable compared to controls", subtitle = "Day 21, 400 ms exposure time, error bar = SEM", fill = "Groups")+
  geom_bar(position = 'dodge', stat = 'summary', fun = 'mean') +
    stat_summary(fun.data = mean_se,  
                 geom = "errorbar", position = 'dodge', colour = "red") + 
  geom_dotplot(data = live_dat2simp, aes(x = Groups.f, y = avg.int),
      alpha = .8, position_dodge(width = 0.9), binaxis = "y", dotsize = 0.6,
      stackdir = "center", stackratio = 1) +
  theme(legend.position="none")

plot2

summary(mod.pff.slope <- lm(avg.int~dose, data = live_dat.21.pff))
plot3 <- ggplot(live_dat.21.pff, aes(x=dose, y=avg.int, fill=dose)) +
  geom_point() + geom_smooth(method = "lm") +
  theme(legend.position="none") +
  ylab("Average Pixel Intensity (PFF tag)") +
  xlab("Dose (µg)") + 
  labs(
    title = "Relationship of PFF dosage to Average Pixel Intensity (Day 21)") +
  theme(text = element_text(size = 14))
  #scale_fill_viridis_d()

plot3

```

```{r percentile comp}
# Plot values
comp1 <- ggplot(live_dat.t.pff, aes(x=day.no,y=max.int,fill=day.no))+
  geom_bar(position = 'dodge', stat = 'summary', fun = 'mean') +
  labs(title = "Maximum pixel intensity",
       x = "Day",
       y = "Maximum Pixel Intensity (ATTO-590)") + 
  geom_dotplot(data = live_dat.t.pff, aes(x = day.no, y = max.int),
               alpha = .8, position_dodge(width = 0.9), binaxis = "y", dotsize = 0.5,
               stackdir = "center", stackratio = 1) +
    geom_signif(comparisons=list(c('7', '9')), annotations="***",
              y_position = 150, tip_length = 0.05, size = 0.8, textsize = 8, vjust=0.4) +
    geom_signif(comparisons=list(c('14','21')), annotations="**",
              y_position = 150, tip_length = 0.05, size = 0.8, textsize = 8, vjust=0.4) +
  stat_summary(fun.data = mean_se,  
    geom = "errorbar", position = 'dodge', colour = "black") +
  theme(legend.position="none",
        text = element_text(size = 14),
        axis.title.x = element_blank()) +
  coord_cartesian(ylim = c(0,165)) +
  scale_fill_viridis_d()

comp1 + labs(title = "Maximum pixel intensity over time (PFF Group)")

comp2 <- ggplot(live_dat.t.pff, aes(x=day.no,y=avg.int,fill=day.no))+
  geom_bar(position = 'dodge', stat = 'summary', fun = 'mean') +
    labs(title = "Average pixel intensity",
       x = "Group",
       y = "Average Pixel Intensity (ATTO-590)") + 
  geom_dotplot(data = live_dat.t.pff, aes(x = day.no, y = avg.int),
               alpha = .8, position_dodge(width = 0.9), binaxis = "y", dotsize = 0.1,
               stackdir = "center", stackratio = 1) +
    geom_signif(comparisons=list(c("PFF", "Control")), annotations="*",
              y_position = 7.8, tip_length = 0.05, vjust=0.4) +
    geom_signif(comparisons=list(c("2XPFF", "Control")), annotations="***",
              y_position = 7.8, tip_length = 0.05, vjust=0.4) +
    geom_signif(comparisons=list(c("4XPFF", "Control")), annotations="***",
              y_position = 7.8, tip_length = 0.05, vjust=0.4) +
  stat_summary(fun.data = mean_se,  
    geom = "errorbar", position = 'dodge', colour = "red")  +
  theme(legend.position="none")


plot_grid(comp2,comp1)


plot4 <- ggplot(live_dat21, aes(x=Groups.f, y=max.int, fill=Groups.f)) +
  ylab("Maximum intensity (unit8)") + xlab("Groups") + 
  labs(title = "Preformed fibril fluorescent tag is measurable compared to controls", subtitle = "Day 21, 400 ms exposure time, error bar = SEM", fill = "Groups")+
  geom_bar(position = 'dodge', stat = 'summary', fun = 'mean') +
    stat_summary(fun.data = mean_se,  
                 geom = "errorbar", position = 'dodge', colour = "red") + 
  geom_dotplot(data = live_dat21, aes(x = Groups.f, y = max.int),
      alpha = .8, position_dodge(width = 0.9), binaxis = "y", dotsize = 0.2,
      stackdir = "center", stackratio = 1) +
  theme(legend.position="none")

plot4

plot5 <- ggplot(live_dat21, aes(x=Groups.f, y=perc90, fill=Groups.f)) +
  ylab("Maximum intensity (unit8)") + xlab("Groups") + 
  labs(title = "Preformed fibril fluorescent tag is measurable compared to controls", subtitle = "Day 21, 400 ms exposure time, error bar = SEM", fill = "Groups")+
  geom_bar(position = 'dodge', stat = 'summary', fun = 'mean') +
    stat_summary(fun.data = mean_se,  
                 geom = "errorbar", position = 'dodge', colour = "red") + 
  geom_dotplot(data = live_dat21, aes(x = Groups.f, y = perc90),
      alpha = .8, position_dodge(width = 0.9), binaxis = "y", dotsize = 0.2,
      stackdir = "center", stackratio = 1) +
  theme(legend.position="none")

plot5

np.pff.t <- pairwise.wilcox.test(live_dat.t.pff$max.int, live_dat.t.pff$day, p.adjust.method = "BH")
np.pff.t
```


```{r timeplots}
# Identify sample size that is analyzed. N = experimental replicate, m = wells, d = number of images
# length(unique(g_cone$rat))

# For each level of information, identify the number of unique values (e.g., attempt number) across each day and experimental group
n.d21 <- aggregate(data = live_dat.t,
                          rat.no ~ day + group,
                          function(rat.no) length(unique(rat.no)))
m.d21 <- aggregate(data = live_dat.t,
                          w.id ~ day + group,
                          function(w.id) length(unique(w.id)))
i.d21  <- aggregate(data = live_dat.t,
                          image.no ~ day + group,
                          function(image.no) length(unique(image.no)))

live.replicates <- cbind(n.d21[,1],n.d21[,3],m.d21[,3],i.d21[,3])
colnames(live.replicates) <- c("Day","N","m","no.images")
rownames(live.replicates) <-m.d21[,2]
live.replicates

time_plot <-ggplot(live_dat.t, aes(x=day, y=avg.int,color=Groups.f)) +
  geom_line(stat = 'summary', fun = 'mean', aes(color=Groups.f)) +
  stat_summary(fun.data = mean_se, geom = "errorbar",aes(color= Groups.f)) +
  labs(title = "Butyrate Modulates PFF Signal Retention over Culture Period", 
       color = "Group") +
  xlab("Time (days)") +
  ylab("Average Pixel Intensity (ATTO-590)") +
  geom_vline(xintercept = 7, colour="gray",linetype = "longdash")

time_plot + 
  annotate("text", x=9, y=0.0045, label= "PFF Application")
```


```{r longitudinal analysis}
# Define variables
live_dat.t$samplesize <- ifelse((live_dat.t$Groups.f == "Control") | (live_dat.t$Groups.f == "But"),1,0)
live_dat.t$dayqual <- ifelse(live_dat.t$day == 7 | live_dat.t$day == 9,1,0)
live_dat.simp <- subset(live_dat.t, samplesize == 0 & dayqual == 0)
live_dat.control <- subset(live_dat.t, samplesize == 1 & dayqual == 0)

# hist(live_dat.simp$avg.int)

# Box-cox transform values
bc.live1 <- boxcox(avg.int ~ Groups.f * day.no, data = live_dat.simp)
(lambda.live1 <- bc.live1$x[which.max(bc.live1$y)])
live_dat.simp$box.int <- ((live_dat.simp$avg.int^lambda.live1-1)/lambda.live1)
skewness(live_dat.simp$box.int, na.rm = FALSE)

# hist(live_dat.simp$box.int)

# Include multiple statistical models
m.live1 <- lm(box.int ~ Groups.f * day.no, data = live_dat.simp)
m.live2 <- lmer(box.int ~ Groups.f * day.no +
                  (1|rat.no), data = live_dat.simp)
m.live3 <- lmer(box.int ~ Groups.f * day.no +
                  (1|rat.no) + (1|rat.no:well), data = live_dat.simp)

# Likelihood ratio tests
lrtest(m.live1,m.live2,m.live3)
lrtest(m.live2,m.live3)

# Re-transform model to allow for back-transformation
boxlive.tran <- make.tran("boxcox", lambda.live1)
m.live.tran <- with(boxlive.tran, 
    lmer(linkfun(box.int) ~ Groups.f * day.no +
         (1|rat.no) + (1|rat.no:well), 
         data = live_dat.simp))

# Present values
pwpm(em.live.tran <- emmeans(m.live.tran, ~ Groups.f * day.no, 
                       type = "response"))
summary(em.live.tran, type = "response")

# Plot values
plot(em.live.tran, comparisons = TRUE, type = "response") + 
  labs(title = "Live Average Pixel Intensity (PFF tag) over time",
       x = "Estimated Marginal Mean (Average Intensity PFF tag)",
       y = "Groups") +
    scale_y_discrete(labels = c(
      "PFF Day 14", "PFF + But Day 14",
      "PFF Day 21", "PFF + But Day 21")) +
  theme(text = element_text(size = 14),
          axis.title.y = element_blank())

```


```{r ranked sum anova}
#Nonparametric
aov.21a <- kruskal.test(avg.int ~ Groups.f, data = live_dat21)
aov.21a
np.d21 <- pairwise.wilcox.test(live_dat21$avg.int, live_dat21$Groups.f,p.adjust.method = "BH")
np.d21

#Box Cox
bc.live.dat21 <- boxcox(live_dat21$avg.int ~ live_dat21$Groups.f)
(lambda <- bc.live.dat21$x[which.max(bc.live.dat21$y)])
live_dat21$box.avg.int <- ((live_dat21$avg.int^lambda-1)/lambda)

#Linear Model
live.21.tran <- make.tran("boxcox", -0.3434343)
m.box.live21 <- with(live.21.tran, 
    lmer(linkfun(avg.int) ~ Groups.f + (1|rat.no), data = live_dat21))
plot(m.box.live21)
qqnorm(resid(m.box.live21))
qqline(resid(m.box.live21))

#Interpretation
em.live.box.d21 <- emmeans(m.box.live21, ~ Groups.f, type = "response")
summary(em.live.box.d21)
pwpm(em.live.box.d21, adjust = "fdr")
```


